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£h ! We report on a series of tests of agreement between three types of N-body 

t— i 1 simulations: PM, P 3 M, and Tree codes. We find good agreement in both the 

individual and the statistical properties only on scales larger than the mean 
interparticle separation. As a result, we question most numerical results at and 
■ below below galaxy scales, either concerning primordial dark matter or baryonic 

matter coupled to it by gravitation. 

, Subject headings: cosmology:miscellaneous-gravitation-hydrodynamics- 

^ e— | methods: numerical-dark matter 

6 ■ 
& : 

c/3 . 1. Introduction 

; If the dark matter in the Universe is some sort of weakly interacting particle, or 

even a population of primordial black holes, then artificial discreteness effects of all kinds 
should be throughly suppressed in numerical simulations of its gravitational clustering. 
Their contribution both to errors in the initial conditions and two-body scattering in the 
dynamical evolution are negligible in the Universe and therefore should be negligible in the 
simulation. If not, the results (including those of baryons, entrained in the gravitational 
field) cannot be relied upon to be consequences of physical initial conditions - regardless of 
whether they resemble astronomical observations. 

It has been argued, based on a variety of numerical studies, that there are real 
problems with the HFLMR (High Force Low Mass Resolution) approach which dominates 



1 Current Address: Hewlett-Packard Company, High Performance Computing Division, 20 Perimeter 
Summit Blvd, MS 1904, Atlanta, GA 30319-1417 

2 Department of Physics and Astronomy, University of Kansas, Lawrence, KS 66045 



- 2 - 



the cosmological N-body methodology. (Peebles, et al. 1989; Melott 1990; Suisalu and 
Saar 1995; Kuhlman et al. 1996; Melott et al. 1997; Park 1997). It is clear from these 
studies that there are unphysical scattering processes taking place in HFLMR codes, but 
the precise effects are not clear. Recently Craig (1997), Klypin et al (1997), and Moore et 
al. (1997) have argued convincingly that the central regions of dark matter halos have been 
seriously misrepresented by numerical results due to lack of resolution. We present here 
some preliminary results of a much larger study (Splinter et al. 1998) which explores the 
limitations of mass resolution in cosmological clustering simulations. 

We use PM (Hockney & Eastwood 1988; Melott 1981, 1986), AP 3 M (Couchman 1991), 
and Tree (Suginohara et al. 1991; Suto 1993) codes. The P 3 M code had adaptive smoothing 
turned off since we wish to compare a standard P 3 M method. The Tree runs use the fixed 
smoothing length in comoving coordinates, and we set a tolerance parameter 9 = 0.2 which 
is considerably smaller (and thus more accurate) than conventional choices (9 = 0.5 ~ 0.75). 
9 merely controls how far the tree expansion is carried, and thus the accuracy of long-range 
forces. 

The initial power spectrum in all cases was P(k) oc k' 1 up to some cutoff, in 
most cases at the Nyquist frequency k = 16kf dictated by the runs with the fewest 
particles. Realization of the corresponding density field was generated using the Zel'dovich 
approximation (Zel'dovich 1970) to perturb the particles from their initial lattice 
(Doroshkevich et al. 1980). All the models are evolved in the Einstein - de Sitter universe 
(Qo — 1)- The comparisons were performed at three different epochs when the nonlinear 
wavenumber k n \ becomes 16k{, 8k{, and 4k{, where k n \ = k n \(A) is defined as 



In the above A denotes the expansion factor (unity at the initial condition), and these 
values of k n \ correspond to the epochs A = 22.36, 42.13, and 92.20, respectively. The latest 
moment we studied corresponds to nonlinearity on the largest scale which does not suffer 
from finite- volume boundary condition problems (Ryden & Gramman 1991; Kauffmann 
& Melott 1992). The specific runs presented here and the model parameters are shown 
in Table 1. We note that PM codes, which have been extensively used in most physical 
applications with large numbers of particles, are much faster than the other two types. 
Thus our 128 3 PM runs took much less CPU time than even the 32 3 Tree or P 3 M runs. 
The typical limitation on PM runs is memory or disk space, while CPU time is the typical 
limitation on P 3 M or Tree runs, a is the absolute scale of force resolution. We define 
e = anl^ so that e = 1 corresponds to smoothing at the mean interparticle separation. 
HFLMR codes usually run with e < 1, but we will examine their behavior over a range in N 
and e, pushing them toward e = 1 by increasing the number of particles while keeping the 
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absolute scale of force resolution a constant. So far as we know, this crucial experiment has 
not previously been done with HFLMR codes. 

Our primary strategy is therefore to highlight the largely unexplored mass resolution 
issue by varying the number of particles while keeping a constant (it is not exactly constant 
because in fact the shape of the short-range softening function is different in all three codes). 
Within a given code a will be constant, so we can spot trends. Between codes softening will 
be of comparable size. We can define r 50 as the radius where the force drops to 50% of the 
Newtonian value. For the P 3 M code, = 0.92a. For the Tree code = 0.87a. For our 
PM code, r 50 = 0.95 grid unit, albeit with considerable scatter in the softened zone. 

In most cases we set a = 1 (in units where the box size is 128) in both the Tree and 
P M codes, and run PM with a 128 3 mesh. (The P 3 M runs had one mesh cell per particle, 
but this is not the factor determining force resolution there.) It is typical of HFLMR codes 
to have a considerably less than l sep , the mean interparticle separation over the entire 
simulation box. In most uses of PM codes, the opposite is true (although in gravitational 
applications it has become customary to have e as small as 1 or 0.5). 

It is important to note that one cannot represent initial power to higher than the 
Nyquist wavenumber of the particles or the mesh of the FFT used to impose the initial 
conditions, whichever is worse (in fact the latter is rarely the problem in cosmology). 
Therefore most of our runs only have initial power up to k c = 16kf, so that comparisons of 
only the effects of divergent numerical integration can be done. We show here results from a 
full 128 3 PM simulation with this power spectrum, and 32 3 simulations with the same force 
resolution, different mass resolution, and the same initial conditions. In order to explore 
the effects of having initial power at higher wave-numbers which is only possible with more 
particles, we have one PM run with k c = 64/cf. 

We now have for the first time a large number of runs using 3 different codes with 
identical initial conditions, identical force resolution, but varying mass resolution. In 
particular the number of particles varies widely enough to study the same physical system 
with values of e typical of HFLMR codes as well as to study the same system with a PM 
code with similar force resolution and matching mass resolution (e = 1). 

It is important to note time-step limitations. Elementary principles of numerical 
stability require that no particle move more than a fraction of a softening length a, or about 
one-half mesh unit for PM, in a single time-step. This condition was amply enforced for 
all three codes. 
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2. Results 

Figure 1 shows power spectra computed on the 128 3 mesh appropriate to the force 
resolution, normalized such that a Poisson distribution with 128 3 particles has P (k) ~ 1. 
The dotted line shows the spectrum of the 32 3 initial conditions. It looks unusual because 
such spectra are not usually shown out to the force resolution scale - even though later 
results usually are shown to this scale. Other types of unperturbed particle arrangements 
than our lattice will have other features, but the overall amplitude is constrained to be 
similar. A full complement of 128 3 particles would allow the power law to continue across 
the plot. 

It should be noted that in all HFLMR codes there are features due to particle 
discreteness present beyond the particle Nyquist frequency, which are resolved and evolved. 
Whatever the particle arrangement, they are not random phase and cannot evolve as the 
initial conditions do since they are tied to the particles. 

Above this we see spectra for three evolved stages. The first (corresponding to the 
onset of nonlinearity at the mass resolution scale of the 32 3 particle simulations) and the 
last (without serious boundary condition problems) of these are shown on the left side as 
ratios to the power in one of them, in order to clarify differences. It can be seen that in 
the first of these that the spectra agree well at low k (linear theory) and very high k (pure 
mode coupling) frequencies, but disagree between. In the range of 10 < k < 20 (in units of 
the fundamental kf), there are three classes of spectra: those with 32 3 particles (bottom 
group), those with the same initial conditions but more particles (middle group) and the 
one which had full initial power down to k = 64 (top group). There is nearly a factor of two 
range here. (In this and following measures we compute spectra from an initially identical 
32 3 subset of the 128 3 run to suppress sampling differences). By the final stage (as evolved 
as is safe with periodic boundary conditions) the differences have moved to high k. If only 
force resolution mattered, these would be identical to k — 64. In fact, they diverge around 
k = 20, close to the limit set by mass resolution. We cannot say which run is correct, and 
we did not even find a monotonic trend with mass resolution. 

An arbitrary density field is described by both the amplitudes and phases of its 
Fourier coefficients. We now examine the phase agreement between various codes. Figure 
2 examines < cos 9 >, where 9 is the difference in phase angle between the same Fourier 
coefficient in different simulations and the average is over all wave-numbers with the same 
amplitude Simulations which agree have < cos 9 >= 1; totally uncorrected results 
have < cos 9 >= 0. The open circles show the 128 3 PM and P 3 M runs with identical initial 
conditions compared. (We do not have a 128 3 Tree run but results with 64 3 and 32 3 runs 
shown in Splinter et al. 1998 support the same trend.) The bottom group (filled squares) 
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shows the 128 3 particle PM run with full initial conditions down to k = 64 compared with 
all other runs. In order to reduce confusion, we do not show our other results, which lie 
in between. Thus (a) As the number of particles are increased, with force resolution and 
initial conditions held constant, different codes agree better, (b) Even when they agree 
(given identical initial conditions), none of them agree strongly on small scales with a run 
that had the full range of proper initial conditions as allowed by having more particles. 

In Table 2, we show the density correlation coefficient k±2 = ^^72 where 5 are the 
density contrasts in two simulations and a their RMS values. Above the diagonal we show 
values on the 128 3 mesh; below the 32 3 mesh. Below the diagonal, results approach k — 1, 
indicating near perfect agreement. Above the diagonal, again recalling that the k c = 16 PM 
run, the Tree, and P 3 M runs have the same initial conditions, we see that these runs agree 
best when the 128 3 runs are compared. Lastly, the k c = 64 PM run does not agree strongly 
with any of the others. Again we conclude that integration errors are greatly reduced when 
more particles are included, but that the inclusion of perturbations on small scales is also 
important. The precise effect of ignoring them is spectrum-dependent. 



3. Does It Matter? 

The reader is referred to the center and bottom center images of Figure 7 in Beacom 
et al. (1991) for another example of the importance of initial power on small scales even 
when it is deep in the nonlinear regime. There, as well as in Melott & Shandarin (1993), we 
emphasized the effectiveness of the transfer of power from long to short waves. The general 
position and orientation of objects is determined by initial perturbations on that comoving 
scale and larger, so smaller perturbations are ignorable for this purpose. See also Little et 
al. (1991), Evrard and Crone (1992). However, here we stress that the internal structure 
of these objects will vary depending on smaller-scale perturbations. If we wish to study 
that internal structure, the smaller-scale initial perturbations must be present and properly 
evolved. 

We have shown that even with nearly identical force resolution, different N-body codes 
differ in their results below the mean interparticle separation. They converge either by 
smoothing on this larger scale or by putting in more particles so that the scales are the 
same. Codes which give different answers cannot all be correct. On the other hand, it might 
be that the results on small scales are statistically equivalent, in the sense that quantities 
of interest computed from the ensembles are the same. This turns out not to be true. Even 
the autocorrelation function is affected by mass resolution; phase sensitive measures are 
more strongly affected. We present here only one simple set of measures, based on the mass 
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density. 

The percolation code of Yess and Shandarin (1996) is used to identify connected 
regions on the 128 3 mesh above a given density threshold. In Figure 3 we show the ratio 
of total mass in regions above a given density contrast to the value found in the k c = 64 
PM simulation at the first and last of our evolved stages. The situation is time-dependent 
and complex to describe. Clearly there are major differences, but the codes agree in low 
density regions. The total mass at the high density limit varies by an order of magnitude. 
At the earlier stage all curves were below unity in the region of 5 ~ 50, suggesting that 
a generation of structures were missed by the absence of correct initial conditions found 
only in one run. Later we find that within a code type, a higher particle density results 
in higher peak densities. Clearly, high density regions are not trustworthy! We found that 
with binning on the 32 3 mesh (not shown) the curves agree well; the maximum difference 
there is about 25%. 

We have examined infall velocities and velocity dispersions (Splinter et al. 1998) and 
found measurable differences. All the differences we find are systematic errors, so no error 
bars are shown. The differences should be measured against the desired accuracy one wants 
to get from the computations. 

We cannot easily compare our results with others, because no study has included our 
variety of codes, mass resolution, and inclusion of normally absent small-scale perturbations. 
The most similar work is the unpublished study coordinated by D.H. Weinberg. An analysis 
similar to our velocity studies (not shown) produced very similar results. Efstathiou et al. 
(1988) studied evolution of power-law initial conditions in a P 3 M code, using primarily 
low-order or averaged statistics. They did find fluctuations of order 50% in the rescaled 
multiplicity function (their Figure 9); while not equivalent, this is roughly compatible with 
the nature of our results in Figure 3. 

We believe the safest course to follow is to restrict attention to scales above the mean 
particle separation. Given current computer technology and volumes large enough to 
respect boundary conditions, this implies scales of order 100 to 200 kpc. With nested-grid 
codes (eg. Splinter 1996) or other schemes to approximate an external zone, this may be 
considerably improved. However, we caution that careful testing is needed in the nonlinear 
regime. It is not enough to produce something that resembles our Universe; we must 
have confidence that it is a consequence of the initial conditions that were supposed to be 
modeled. 

At the present time, we are quite skeptical of nearly all numerical results on early 
galaxy formation and the inner parts of dark matter halos, as they are typically below the 
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mean interparticle separation of the simulations in question. 
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Table 1. Model Parameters for the Test Cases 



Code N 



PM k c = 16 128 3 1.0 1.0 

k c = 64 12S 3 1.0 1.0 

P 3 M k c = 16 32 3 4.0 0.25 

k c = 16 12S 3 1.0 1.0 

Tree k c = 16 32 3 4.0 0.25 



a ^sep i s the mean particle separation in grid cell units. 

b ey 01 . ce is in units of the mean particle separation, T sep . Note that for all runs here a = £f orce X l B 



Table 2. Cross-Correlations at Final Stage 





fc c = 16 


k c = 64 


128 3 e = 1.0 


e = 0.25 


e = 0.25 


PM (k c = 16) 




0.65 


0.83 


0.64 


0.52 


PM (k c = 64) 


0.96 




0.63 


0.57 


0.47 


P 3 M 128 3 (e = 1.0) 


0.99 


0.95 




0.66 


0.50 


P 3 M(e = 0.25) 


0.96 


0.94 


0.97 




0.70 


Trec(e = 0.25) 


0.92 


0.90 


0.92 


0.96 





-10- 




Fig. 1. — Right panel: The power spectrum of our initial conditions and three successive 
evolved stages for all our simulations, evaluated from 32 3 particles on a 128 3 mesh (their force 
resolution). The normalization is such that a Poisson distribution of 128 3 particles would 
converge to P = 1. The light solid line is the k c = 16A; f PM run; the heavy solid line is the 
k c = 64fcf PM run. Other heavy lines are P 3 M runs and light lines are tree code runs. The 
dotted lines are the initial conditions for the 32 3 particle runs; the spikes are discreteness 
effects due to the finite number of particles. Other lines: Longdash is the P 3 M run with 
128 3 particles, and the longdash-dot lines are the P 3 M and Tree runs with 32 3 particles. Left 
panel: The ratio of the power in a given model to that in our fiducial k c = 64/cf PM run at 
the first (bottom) and last (top) evolved stage. 
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Fig. 2. — The successive plots contain all the data for the averaged phase agreement between 
all of the simulations runs at the same stage; < cos 9 > is defined in the text and is 1 
for agreement and for uncorrelated phases. The filled squares indicate cross-correlation of 
other runs against the PM run which continued the power-law perturbations to wave-numbers 
impossible for small numbers of particles (k c = 6Ak{). Open circles represent cross-correlation 
between the 128 3 P 3 M run and the PM run with k c = 16k{ (also 128 3 particles). All other 
possible cross-correlations here (not shown) lie between these extremes. As before, the top 
panel is the last evolved stage. 
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Fig. 3. — The lines (same types as in Figure 1) show the amount of mass in regions of 
densities greater than the density threshold shown on the horizontal. The densities were 
calculated on the 128 3 mesh. The values are shown as ratios to the fiducial PM model with 
N = 128 3 , and k c = 64, at our earliest evolved stage (bottom) and our latest (top). 



